Industrial Production Index
Qian Yi
4/12/2020
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(fGarch) ## Loading required package: timeDate## Loading required package: timeSeries## Loading required package: fBasicslibrary(dynlm)## Loading required package: zoo##
## Attaching package: 'zoo'## The following object is masked from 'package:timeSeries':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(tseries)## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(astsa) ##
## Attaching package: 'astsa'## The following object is masked from 'package:fBasics':
##
## nyselibrary(xts)##
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, lastlibrary(fpp)## Loading required package: forecast## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.fracdiff fracdiff
## residuals.fracdiff fracdiff##
## Attaching package: 'forecast'## The following object is masked from 'package:astsa':
##
## gas## Loading required package: fma##
## Attaching package: 'fma'## The following objects are masked from 'package:astsa':
##
## chicken, sales## Loading required package: expsmooth## Loading required package: lmtest##
## Attaching package: 'fpp'## The following object is masked from 'package:astsa':
##
## oillibrary(TSA)## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(forecast)
library(dplR)
library(tseries)Read the dataset
ifs <- read.csv("IFS.csv",stringsAsFactors =FALSE)
colnames(ifs)=c("Country_Name", "Country_Code", "Indicator_Name", "Indicator_Code", "Time_Period", "Value", "Status")
# remove Status column
ifs=ifs[, c(1:6)]
# take a look at the dataset
head(ifs)## Country_Name Country_Code
## 1 United States 111
## 2 United States 111
## 3 United States 111
## 4 United States 111
## 5 United States 111
## 6 United States 111
## Indicator_Name
## 1 International Reserves, Official Reserve Assets, US Dollars
## 2 International Reserves, Official Reserve Assets, US Dollars
## 3 International Reserves, Official Reserve Assets, US Dollars
## 4 International Reserves, Official Reserve Assets, US Dollars
## 5 International Reserves, Official Reserve Assets, US Dollars
## 6 International Reserves, Official Reserve Assets, US Dollars
## Indicator_Code Time_Period Value
## 1 RAFA_USD 1980M5 21917557876
## 2 RAFA_USD 1980M9 22994402034
## 3 RAFA_USD 1980M11 25673025848
## 4 RAFA_USD 1981M4 29692779278
## 5 RAFA_USD 1981M9 29715202935
## 6 RAFA_USD 1982M4 31555856317# dim(ifs)
print(colnames(ifs))## [1] "Country_Name" "Country_Code" "Indicator_Name" "Indicator_Code"
## [5] "Time_Period" "Value"United　States
Research Indicator: Economic Activity, Industrial Production, Index
Indicator Code: AIP_IX
Country: USA
df1_us=subset(ifs, (Indicator_Code=='AIP_IX') & (Country_Name=='United States'))now the time didn’t apear as ascending order
df1_us$Time_Period <- gsub("M", "/", df1_us$Time_Period)
df1_us['date'] = '/01'
df1_us$Time_Period <- paste(df1_us$Time_Period, df1_us$date, sep="")
df1_us$Time_Period = as.Date(df1_us$Time_Period, "%Y/%m/%d")
df1_us = df1_us[order(df1_us$Time_Period),]
df1_us = subset(df1_us, select = -c(date))
rownames(df1_us)=NULLOnly keep value and Time_Period
# Research Indicator: Economic Activity, Industrial Production, Index
# Indicator Code: AIP_IX
# Country: USA
# the original data
df11=ts(df1_us$Value, frequency = 12, start = c(1980, 1))
plot(df11, main="United States: Industrial Production Index", xlab='year', ylab='US IPI')
(1) decompose Analysis
plot(decompose(df11))
(2) Spectral Analysis
sp=spectrum(df11, log='no') 
# find the cycle with predominant frequecny
sp$fre[which.max(sp$spec)] # 0.025## [1] 0.025# find the cycle of predominant frequecny
1/sp$fre[which.max(sp$spec)] # 40 month## [1] 40max(sp$spec) #estimate of spectral density at predominant period which is close to 326.3907.## [1] 326.3907U = qchisq(.025,2)
L = qchisq(.975,2)
2*max(sp$spec)/L #lower bound for spectral density## [1] 88.479622*max(sp$spec)/U #upper bound for spectral density## [1] 12891.74# The 95% confidence interval is [88.47962, 12891.74](3) ARIMA model
fit1 <- auto.arima(df11, trace = TRUE) # Best model: ARIMA(4,0,1)(2,1,1)[12] with drift # AIC: 1022.693##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2)(1,1,1)[12] with drift : 987.9367
## ARIMA(0,0,0)(0,1,0)[12] with drift : 2487.57
## ARIMA(1,0,0)(1,1,0)[12] with drift : 1102.752
## ARIMA(0,0,1)(0,1,1)[12] with drift : 1972.828
## ARIMA(0,0,0)(0,1,0)[12] : 2566.802
## ARIMA(2,0,2)(0,1,1)[12] with drift : 999.8981
## ARIMA(2,0,2)(1,1,0)[12] with drift : 1029.682
## ARIMA(2,0,2)(2,1,1)[12] with drift : 985.0017
## ARIMA(2,0,2)(2,1,0)[12] with drift : 989.1653
## ARIMA(2,0,2)(2,1,2)[12] with drift : 979.5924
## ARIMA(2,0,2)(1,1,2)[12] with drift : 980.8966
## ARIMA(1,0,2)(2,1,2)[12] with drift : 1022.589
## ARIMA(2,0,1)(2,1,2)[12] with drift : 1003.046
## ARIMA(3,0,2)(2,1,2)[12] with drift : 1041.516
## ARIMA(2,0,3)(2,1,2)[12] with drift : 978.6119
## ARIMA(2,0,3)(1,1,2)[12] with drift : 979.8452
## ARIMA(2,0,3)(2,1,1)[12] with drift : 983.9795
## ARIMA(2,0,3)(1,1,1)[12] with drift : 987.6149
## ARIMA(1,0,3)(2,1,2)[12] with drift : 999.6705
## ARIMA(3,0,3)(2,1,2)[12] with drift : 971.2854
## ARIMA(3,0,3)(1,1,2)[12] with drift : 973.9898
## ARIMA(3,0,3)(2,1,1)[12] with drift : 969.51
## ARIMA(3,0,3)(1,1,1)[12] with drift : 975.4823
## ARIMA(3,0,3)(2,1,0)[12] with drift : 979.8806
## ARIMA(3,0,3)(1,1,0)[12] with drift : 1096.151
## ARIMA(3,0,2)(2,1,1)[12] with drift : 1044.916
## ARIMA(4,0,3)(2,1,1)[12] with drift : 969.8308
## ARIMA(3,0,4)(2,1,1)[12] with drift : 971.0494
## ARIMA(2,0,4)(2,1,1)[12] with drift : 985.9654
## ARIMA(4,0,2)(2,1,1)[12] with drift : 967.7735
## ARIMA(4,0,2)(1,1,1)[12] with drift : 976.2705
## ARIMA(4,0,2)(2,1,0)[12] with drift : 976.1243
## ARIMA(4,0,2)(2,1,2)[12] with drift : 969.0887
## ARIMA(4,0,2)(1,1,0)[12] with drift : 1022.999
## ARIMA(4,0,2)(1,1,2)[12] with drift : 975.0806
## ARIMA(4,0,1)(2,1,1)[12] with drift : 966.197
## ARIMA(4,0,1)(1,1,1)[12] with drift : 974.9985
## ARIMA(4,0,1)(2,1,0)[12] with drift : 974.9076
## ARIMA(4,0,1)(2,1,2)[12] with drift : 967.5903
## ARIMA(4,0,1)(1,1,0)[12] with drift : 1021.926
## ARIMA(4,0,1)(1,1,2)[12] with drift : 973.5646
## ARIMA(3,0,1)(2,1,1)[12] with drift : 974.9792
## ARIMA(4,0,0)(2,1,1)[12] with drift : 974.2169
## ARIMA(5,0,1)(2,1,1)[12] with drift : 968.7653
## ARIMA(3,0,0)(2,1,1)[12] with drift : 1016.15
## ARIMA(5,0,0)(2,1,1)[12] with drift : 971.5388
## ARIMA(5,0,2)(2,1,1)[12] with drift : 970.6539
## ARIMA(4,0,1)(2,1,1)[12] : 970.5866
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(4,0,1)(2,1,1)[12] with drift : 1022.693
##
## Best model: ARIMA(4,0,1)(2,1,1)[12] with driftevaluate the model
sarima(df11, 4,0,1,2,1,1,12) # accoring to the disgnose plots, the model is a good fit ## initial value 1.283466
## iter 2 value 0.336696
## iter 3 value 0.234533
## iter 4 value 0.062828
## iter 5 value 0.004002
## iter 6 value -0.085041
## iter 7 value -0.129895
## iter 8 value -0.218406
## iter 9 value -0.285143
## iter 10 value -0.299936
## iter 11 value -0.302313
## iter 12 value -0.305571
## iter 13 value -0.320862
## iter 14 value -0.328682
## iter 15 value -0.340570
## iter 16 value -0.344303
## iter 17 value -0.348020
## iter 18 value -0.352979
## iter 19 value -0.361265
## iter 20 value -0.365733
## iter 21 value -0.365747
## iter 22 value -0.366222
## iter 23 value -0.366602
## iter 24 value -0.366917
## iter 25 value -0.367217
## iter 26 value -0.367664
## iter 27 value -0.368020
## iter 28 value -0.368310
## iter 29 value -0.368590
## iter 30 value -0.368807
## iter 31 value -0.368887
## iter 32 value -0.368905
## iter 33 value -0.368920
## iter 34 value -0.368947
## iter 35 value -0.368975
## iter 36 value -0.368989
## iter 37 value -0.368991
## iter 38 value -0.368991
## iter 39 value -0.368991
## iter 39 value -0.368991
## final value -0.368991
## converged
## initial value -0.342462
## iter 2 value -0.342870
## iter 3 value -0.343386
## iter 4 value -0.343408
## iter 5 value -0.343445
## iter 6 value -0.343489
## iter 7 value -0.343524
## iter 8 value -0.343534
## iter 9 value -0.343543
## iter 10 value -0.343559
## iter 11 value -0.343579
## iter 12 value -0.343592
## iter 13 value -0.343596
## iter 14 value -0.343598
## iter 15 value -0.343601
## iter 16 value -0.343603
## iter 17 value -0.343604
## iter 18 value -0.343604
## iter 19 value -0.343605
## iter 20 value -0.343605
## iter 21 value -0.343605
## iter 22 value -0.343605
## iter 23 value -0.343605
## iter 24 value -0.343606
## iter 25 value -0.343606
## iter 26 value -0.343606
## iter 27 value -0.343606
## iter 28 value -0.343606
## iter 29 value -0.343606
## iter 30 value -0.343606
## iter 30 value -0.343606
## iter 30 value -0.343606
## final value -0.343606
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 sar1 sar2 sma1
## 1.3971 -0.2448 0.0066 -0.1750 -0.4334 -0.1474 -0.1339 -0.4812
## s.e. 0.1644 0.1835 0.0907 0.0701 0.1669 0.1127 0.0767 0.1086
## constant
## 0.1231
## s.e. 0.0390
##
## sigma^2 estimated as 0.4928: log likelihood = -501.1, aic = 1022.21
##
## $degrees_of_freedom
## [1] 457
##
## $ttable
## Estimate SE t.value p.value
## ar1 1.3971 0.1644 8.4990 0.0000
## ar2 -0.2448 0.1835 -1.3338 0.1829
## ar3 0.0066 0.0907 0.0726 0.9421
## ar4 -0.1750 0.0701 -2.4959 0.0129
## ma1 -0.4334 0.1669 -2.5967 0.0097
## sar1 -0.1474 0.1127 -1.3076 0.1917
## sar2 -0.1339 0.0767 -1.7467 0.0814
## sma1 -0.4812 0.1086 -4.4301 0.0000
## constant 0.1231 0.0390 3.1602 0.0017
##
## $AIC
## [1] 2.142997
##
## $AICc
## [1] 2.143805
##
## $BIC
## [1] 2.229878fit ARIMA by hand
plot(diff(df11)) # d=1
acf(diff(df11), lag.max=200) # MA:q= 0
pacf(diff(df11), lag.max=200) # AR:p= 4,5
# try ARIMA(4,0,1)
fit1.1 <- Arima(df11, order=c(4, 1, 1), seasonal=c(0,1,0), method="ML",include.drift = TRUE)## Warning in Arima(df11, order = c(4, 1, 1), seasonal = c(0, 1, 0), method =
## "ML", : No drift term fitted as the order of difference is 2 or more.summary(fit1.1)## Series: df11
## ARIMA(4,1,1)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1
## 0.0380 0.1928 0.2395 0.0621 -0.0319
## s.e. 0.5318 0.0466 0.1189 0.1400 0.5322
##
## sigma^2 estimated as 0.7274: log likelihood=-583.42
## AIC=1178.83 AICc=1179.02 BIC=1203.68
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -0.00228956 0.8366436 0.6124725 0.004780768 0.7093855
## MASE ACF1
## Training set 0.1979926 0.00137799checkresiduals(fit1.1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,1,1)(0,1,0)[12]
## Q* = 112.44, df = 19, p-value = 2.776e-15
##
## Model df: 5. Total lags used: 24# AIC=1178.83 AICc=1179.02 BIC=1203.68
fit1.2 <- Arima(df11, order=c(5, 0, 1), seasonal=c(0,1,0), method="ML",include.drift = TRUE)
summary(fit1.2)## Series: df11
## ARIMA(5,0,1)(0,1,0)[12] with drift
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 drift
## 0.7600 0.3868 0.1127 -0.1658 -0.1681 0.1897 0.1215
## s.e. 1.5732 1.5818 0.2922 0.0942 0.3065 1.6364 0.0501
##
## sigma^2 estimated as 0.6837: log likelihood=-570.79
## AIC=1157.57 AICc=1157.89 BIC=1190.73
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 0.005559697 0.8102737 0.6024348 0.009176348 0.6974393
## MASE ACF1
## Training set 0.1947477 -0.01055751checkresiduals(fit1.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(5,0,1)(0,1,0)[12] with drift
## Q* = 102.35, df = 17, p-value = 3.253e-14
##
## Model df: 7. Total lags used: 24# AIC=1157.57 AICc=1157.89 BIC=1190.73
# comparing the by hand models with ARIMA(4,0,1)(2,1,1)[12], the ARIMA(4,0,1)(2,1,1)[12] fits betterUse ARIMA model to predict future index
pred1=forecast(fit1,h=48)
pred1 # The last 2 columns are 95% prediction intervals lower bound and upper bound ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Nov 2019 114.9564 114.0478 115.8649 113.5669 116.3459
## Dec 2019 114.9664 113.7047 116.2282 113.0367 116.8961
## Jan 2020 114.8472 113.2367 116.4577 112.3841 117.3103
## Feb 2020 114.4272 112.4248 116.4297 111.3647 117.4897
## Mar 2020 114.8748 112.5066 117.2430 111.2530 118.4967
## Apr 2020 113.1381 110.4227 115.8535 108.9853 117.2909
## May 2020 113.5341 110.4884 116.5797 108.8762 118.1920
## Jun 2020 116.6900 113.3371 120.0428 111.5623 121.8177
## Jul 2020 114.8741 111.2374 118.5108 109.3122 120.4360
## Aug 2020 117.6623 113.7649 121.5597 111.7018 123.6229
## Sep 2020 116.0626 111.9277 120.1974 109.7388 122.3863
## Oct 2020 115.3294 110.9796 119.6792 108.6769 121.9818
## Nov 2020 115.0089 110.3571 119.6607 107.8945 122.1232
## Dec 2020 115.2870 110.3707 120.2033 107.7682 122.8058
## Jan 2021 115.3925 110.2268 120.5583 107.4922 123.2928
## Feb 2021 115.2023 109.7942 120.6105 106.9313 123.4734
## Mar 2021 115.8583 110.2268 121.4899 107.2456 124.4711
## Apr 2021 114.4128 108.5759 120.2497 105.4861 123.3395
## May 2021 114.8070 108.7824 120.8316 105.5931 124.0209
## Jun 2021 118.1339 111.9396 124.3281 108.6606 127.6072
## Jul 2021 116.4791 110.1327 122.8255 106.7731 126.1850
## Aug 2021 119.4140 112.9320 125.8959 109.5007 129.3273
## Sep 2021 117.9293 111.3274 124.5311 107.8325 128.0260
## Oct 2021 117.2416 110.5343 123.9489 106.9836 127.4996
## Nov 2021 117.0175 110.1625 123.8725 106.5337 127.5013
## Dec 2021 117.2923 110.3086 124.2761 106.6116 127.9731
## Jan 2022 117.4746 110.3696 124.5796 106.6084 128.3408
## Feb 2022 117.2795 110.0563 124.5027 106.2325 128.3264
## Mar 2022 117.9327 110.6005 125.2649 106.7190 129.1463
## Apr 2022 116.3650 108.9327 123.7974 104.9982 127.7318
## May 2022 116.7950 109.2710 124.3190 105.2881 128.3019
## Jun 2022 120.0865 112.4798 127.6931 108.4531 131.7198
## Jul 2022 118.3912 110.7106 126.0718 106.6448 130.1377
## Aug 2022 121.3911 113.6448 129.1375 109.5441 133.2381
## Sep 2022 119.8427 112.0384 127.6470 107.9071 131.7783
## Oct 2022 119.0153 111.1603 126.8703 107.0022 131.0284
## Nov 2022 118.7702 110.8218 126.7185 106.6142 130.9261
## Dec 2022 119.0077 110.9775 127.0379 106.7265 131.2888
## Jan 2023 119.1454 111.0346 127.2561 106.7410 131.5497
## Feb 2023 118.9160 110.7219 127.1101 106.3843 131.4477
## Mar 2023 119.5366 111.2631 127.8101 106.8834 132.1899
## Apr 2023 117.9422 109.5932 126.2912 105.1735 130.7109
## May 2023 118.3609 109.9406 126.7812 105.4832 131.2386
## Jun 2023 121.6281 113.1417 130.1144 108.6493 134.6068
## Jul 2023 119.9103 111.3634 128.4573 106.8389 132.9818
## Aug 2023 122.8740 114.2719 131.4761 109.7182 136.0297
## Sep 2023 121.3124 112.6607 129.9642 108.0808 134.5441
## Oct 2023 120.4925 111.7965 129.1884 107.1931 133.7918autoplot(pred1) #plot with confidence interval 
(4) ARCH/GARCH model
# fit ARIMA model without seasonal term
fit2=auto.arima(df11, seasonal = FALSE, trace = TRUE) # ARIMA(2,1,2)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : 1670.596
## ARIMA(0,1,0) with drift : 1849.828
## ARIMA(1,1,0) with drift : 1745.937
## ARIMA(0,1,1) with drift : 1749.283
## ARIMA(0,1,0) : 1850.404
## ARIMA(1,1,2) with drift : 1734.904
## ARIMA(2,1,1) with drift : 1731.286
## ARIMA(3,1,2) with drift : 1727.289
## ARIMA(2,1,3) with drift : 1634.972
## ARIMA(1,1,3) with drift : 1733.839
## ARIMA(3,1,3) with drift : Inf
## ARIMA(2,1,4) with drift : Inf
## ARIMA(1,1,4) with drift : 1731.388
## ARIMA(3,1,4) with drift : 1637.621
## ARIMA(2,1,3) : 1641.337
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,1,3) with drift : Inf
## ARIMA(3,1,4) with drift : Inf
## ARIMA(2,1,3) : Inf
## ARIMA(2,1,2) with drift : 1670.536
##
## Best model: ARIMA(2,1,2) with driftdf11.diff=diff(df11)
plot(df11.diff, main=" diff data") 
df11.log=log(df11)
plot(df11.log, main=" log data")
di11.difflog=diff(log(df11))
plot(di11.difflog, main=" difflog data") 
acf(di11.difflog, lag.max = 100) # from pacf plot, we could observe the seasonality
pacf(di11.difflog, lag.max = 100)# from pacf plot, AR(1) or AR(2)
acf2(df11^2) # accoring to pacf plot, ARCH(1) or ARCH(2) model 
## ACF PACF
## [1,] 0.99 0.99
## [2,] 0.99 0.21
## [3,] 0.98 -0.03
## [4,] 0.97 -0.08
## [5,] 0.97 0.02
## [6,] 0.96 0.05
## [7,] 0.95 -0.07
## [8,] 0.95 -0.06
## [9,] 0.94 0.08
## [10,] 0.93 -0.02
## [11,] 0.93 -0.06
## [12,] 0.92 0.20
## [13,] 0.91 -0.35
## [14,] 0.90 0.10
## [15,] 0.90 0.01
## [16,] 0.89 -0.01
## [17,] 0.88 0.07
## [18,] 0.88 0.03
## [19,] 0.87 -0.03
## [20,] 0.86 0.04
## [21,] 0.86 0.02
## [22,] 0.85 0.02
## [23,] 0.84 -0.04
## [24,] 0.84 0.09
## [25,] 0.83 -0.22
## [26,] 0.83 0.09
## [27,] 0.82 0.01
## [28,] 0.81 -0.04
## [29,] 0.81 0.03
## [30,] 0.80 -0.01
## [31,] 0.79 0.02
## [32,] 0.79 0.01
## [33,] 0.78 0.01
## [34,] 0.78 0.01
## [35,] 0.77 -0.01
## [36,] 0.77 0.06
## [37,] 0.76 -0.17
## [38,] 0.76 0.02
## [39,] 0.75 0.02
## [40,] 0.74 0.00
## [41,] 0.74 0.03
## [42,] 0.73 0.02
## [43,] 0.73 -0.01
## [44,] 0.72 0.04
## [45,] 0.72 -0.02
## [46,] 0.71 0.00
## [47,] 0.71 0.01
## [48,] 0.70 0.01# fit ARIMA(2,1,2) model, then plot the residuals ACF and PACF
arima212=arima(df11, order=c(2,1,2))
residual_arima212=arima212$residuals
acf2(residual_arima212^2)
## ACF PACF
## [1,] 0.03 0.03
## [2,] 0.12 0.12
## [3,] 0.12 0.11
## [4,] -0.07 -0.09
## [5,] 0.05 0.02
## [6,] -0.10 -0.10
## [7,] 0.08 0.10
## [8,] -0.09 -0.09
## [9,] 0.00 0.02
## [10,] 0.10 0.09
## [11,] 0.08 0.12
## [12,] 0.37 0.34
## [13,] 0.01 -0.02
## [14,] 0.06 -0.06
## [15,] 0.04 -0.04
## [16,] -0.08 -0.04
## [17,] 0.03 0.01
## [18,] -0.08 -0.02
## [19,] 0.02 0.00
## [20,] -0.07 -0.02
## [21,] 0.05 0.08
## [22,] 0.11 0.05
## [23,] 0.04 -0.01
## [24,] 0.29 0.14
## [25,] 0.04 0.04
## [26,] 0.07 0.03
## [27,] 0.10 0.08
## [28,] -0.10 -0.08
## [29,] 0.06 0.03
## [30,] -0.06 0.00
## [31,] 0.02 0.01
## [32,] -0.11 -0.11
## [33,] 0.06 0.05
## [34,] 0.03 -0.05
## [35,] 0.00 0.01
## [36,] 0.32 0.18
## [37,] 0.04 0.04
## [38,] 0.11 0.04
## [39,] 0.09 0.02
## [40,] -0.07 -0.04
## [41,] 0.02 -0.03
## [42,] -0.10 -0.04
## [43,] 0.02 0.01
## [44,] -0.10 -0.01
## [45,] -0.01 -0.01
## [46,] 0.09 0.06
## [47,] 0.01 0.00
## [48,] 0.28 0.07par(mfrow=c(2,2))
acf(residual_arima212, main="ACF of residuals (US)", lag.max=100)
pacf(residual_arima212, main="PACF of residuals (US)", lag.max=100)
acf(residual_arima212^2, main="ACF of squared residuals (US)", lag.max=100)
pacf(residual_arima212^2, main="PACF of squared residuals (US)", lag.max=100)
garch_us1=garch(residual_arima212,order=c(3,3),trace=F)
summary(garch_us1)##
## Call:
## garch(x = residual_arima212, order = c(3, 3), trace = F)
##
## Model:
## GARCH(3,3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.72215 -0.57499 0.05402 0.77642 2.72311
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 1.455e+00 2.754e+00 0.529 0.597
## a1 1.054e-02 6.284e-02 0.168 0.867
## a2 1.015e-01 8.114e-02 1.251 0.211
## a3 6.450e-02 2.579e-01 0.250 0.803
## b1 1.073e-02 2.530e+00 0.004 0.997
## b2 1.471e-15 1.614e+00 0.000 1.000
## b3 3.760e-02 8.583e-01 0.044 0.965
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 7.9274, df = 2, p-value = 0.01899
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.055499, df = 1, p-value = 0.8138AIC(garch_us1)## [1] 1659.3garch_us2=garch(residual_arima212,order=c(0,3),trace=F) # (q,p)
summary(garch_us2)##
## Call:
## garch(x = residual_arima212, order = c(0, 3), trace = F)
##
## Model:
## GARCH(0,3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.66389 -0.57327 0.05375 0.77225 2.69503
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 1.564374 0.183095 8.544 <2e-16 ***
## a1 0.009081 0.060293 0.151 0.880
## a2 0.095794 0.076000 1.260 0.208
## a3 0.070207 0.051845 1.354 0.176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 7.0134, df = 2, p-value = 0.03
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.073098, df = 1, p-value = 0.7869AIC(garch_us2) # the AIC is smaller, so ARCH(3) model is better ## [1] 1652.923United　Kingdom
Research Indicator: Economic Activity, Industrial Production, Index
Indicator Code: AIP_IX
Country: United Kingdom
df1_uk=subset(ifs, (Indicator_Code=='AIP_IX') & (Country_Name=='United Kingdom'))
head(df1_uk)## Country_Name Country_Code
## 103049 United Kingdom 112
## 103050 United Kingdom 112
## 103051 United Kingdom 112
## 103052 United Kingdom 112
## 103053 United Kingdom 112
## 103054 United Kingdom 112
## Indicator_Name Indicator_Code
## 103049 Economic Activity, Industrial Production, Index AIP_IX
## 103050 Economic Activity, Industrial Production, Index AIP_IX
## 103051 Economic Activity, Industrial Production, Index AIP_IX
## 103052 Economic Activity, Industrial Production, Index AIP_IX
## 103053 Economic Activity, Industrial Production, Index AIP_IX
## 103054 Economic Activity, Industrial Production, Index AIP_IX
## Time_Period Value
## 103049 1997M8 96.15385
## 103050 1998M1 106.01578
## 103051 1998M7 108.48126
## 103052 1998M10 113.70809
## 103053 1998M12 108.48126
## 103054 1999M3 120.71006dim(df1_uk)## [1] 474 6change the format of date
df1_uk$Time_Period <- gsub("M", "/", df1_uk$Time_Period)
df1_uk['date'] = '/01'
df1_uk$Time_Period <- paste(df1_uk$Time_Period, df1_uk$date, sep="")
df1_uk$Time_Period = as.Date(df1_uk$Time_Period, "%Y/%m/%d")
df1_uk = df1_uk[order(df1_uk$Time_Period),]
df1_uk = subset(df1_uk, select = -c(date))
rownames(df1_uk)=NULLOnly keep value and Time_Period
# Research Indicator: Economic Activity, Industrial Production, Index
# Indicator Code: AIP_IX
# Country: USA
# the original data
df22=ts(df1_uk$Value, frequency = 12, start = c(1980, 1))
plot(df22, main="United Kingdom: Industrial Production Index")
(1) decompose Analysis
plot(decompose(df22))
(2) Spectral Analysis
sp=spectrum(df22, log='no') 
# find the cycle with predominant frequecny
sp$fre[which.max(sp$spec)] # 0.025## [1] 0.025# find the cycle of predominant frequecny
1/sp$fre[which.max(sp$spec)] # 40 month## [1] 40max(sp$spec) #estimate of spectral density at predominant period which is close to 841.0988.## [1] 841.0988U = qchisq(.025,2)
L = qchisq(.975,2)
2*max(sp$spec)/L #lower bound for spectral density## [1] 228.00932*max(sp$spec)/U #upper bound for spectral density## [1] 33221.63# The 95% confidence interval is [228.0093, 33221.63](3) ARIMA model
fit22 <- auto.arima(df22, trace = TRUE) # Best model: ARIMA(5,1,1)(2,1,2)[12] ##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2)(1,1,1)[12] : 2059.437
## ARIMA(0,1,0)(0,1,0)[12] : 2571.477
## ARIMA(1,1,0)(1,1,0)[12] : 2401.882
## ARIMA(0,1,1)(0,1,1)[12] : 2245.311
## ARIMA(2,1,2)(0,1,1)[12] : 2082.606
## ARIMA(2,1,2)(1,1,0)[12] : 2157.584
## ARIMA(2,1,2)(2,1,1)[12] : 2054.25
## ARIMA(2,1,2)(2,1,0)[12] : 2125.736
## ARIMA(2,1,2)(2,1,2)[12] : 2046.093
## ARIMA(2,1,2)(1,1,2)[12] : 2052.854
## ARIMA(1,1,2)(2,1,2)[12] : 2072.008
## ARIMA(2,1,1)(2,1,2)[12] : 2044.492
## ARIMA(2,1,1)(1,1,2)[12] : 2050.995
## ARIMA(2,1,1)(2,1,1)[12] : 2052.372
## ARIMA(2,1,1)(1,1,1)[12] : 2057.492
## ARIMA(1,1,1)(2,1,2)[12] : 2082.583
## ARIMA(2,1,0)(2,1,2)[12] : 2044.521
## ARIMA(3,1,1)(2,1,2)[12] : 2042.901
## ARIMA(3,1,1)(1,1,2)[12] : 2054.019
## ARIMA(3,1,1)(2,1,1)[12] : 2056.749
## ARIMA(3,1,1)(1,1,1)[12] : 2059.586
## ARIMA(3,1,0)(2,1,2)[12] : 2043.263
## ARIMA(4,1,1)(2,1,2)[12] : 2047.265
## ARIMA(3,1,2)(2,1,2)[12] : 2039.889
## ARIMA(3,1,2)(1,1,2)[12] : 2055.05
## ARIMA(3,1,2)(2,1,1)[12] : 2054.321
## ARIMA(3,1,2)(1,1,1)[12] : 2060.681
## ARIMA(4,1,2)(2,1,2)[12] : 2034.072
## ARIMA(4,1,2)(1,1,2)[12] : Inf
## ARIMA(4,1,2)(2,1,1)[12] : 2055.666
## ARIMA(4,1,2)(1,1,1)[12] : 2061.153
## ARIMA(5,1,2)(2,1,2)[12] : 2034.261
## ARIMA(4,1,3)(2,1,2)[12] : 2036.166
## ARIMA(3,1,3)(2,1,2)[12] : Inf
## ARIMA(5,1,1)(2,1,2)[12] : 2032.826
## ARIMA(5,1,1)(1,1,2)[12] : 2050.585
## ARIMA(5,1,1)(2,1,1)[12] : 2056.553
## ARIMA(5,1,1)(1,1,1)[12] : 2055.263
## ARIMA(5,1,0)(2,1,2)[12] : 2046.034
## ARIMA(4,1,0)(2,1,2)[12] : 2045.187
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(5,1,1)(2,1,2)[12] : 2065.401
##
## Best model: ARIMA(5,1,1)(2,1,2)[12]evaluate the model
sarima(df22, 5,1,1,2,1,2,12) ## initial value 1.434990
## iter 2 value 1.230325
## iter 3 value 0.933385
## iter 4 value 0.874606
## iter 5 value 0.851026
## iter 6 value 0.848742
## iter 7 value 0.845081
## iter 8 value 0.840059
## iter 9 value 0.837000
## iter 10 value 0.836880
## iter 11 value 0.836817
## iter 12 value 0.836764
## iter 13 value 0.836583
## iter 14 value 0.836095
## iter 15 value 0.833691
## iter 16 value 0.833226
## iter 17 value 0.832536
## iter 18 value 0.827505
## iter 19 value 0.822364
## iter 20 value 0.820549
## iter 21 value 0.818534
## iter 22 value 0.818008
## iter 23 value 0.817984
## iter 24 value 0.817974
## iter 25 value 0.817901
## iter 26 value 0.817868
## iter 27 value 0.817826
## iter 28 value 0.817781
## iter 29 value 0.817756
## iter 30 value 0.817708
## iter 31 value 0.817648
## iter 32 value 0.817429
## iter 33 value 0.815713
## iter 34 value 0.814529
## iter 35 value 0.812862
## iter 36 value 0.812569
## iter 37 value 0.811984
## iter 38 value 0.811157
## iter 39 value 0.809988
## iter 40 value 0.808249
## iter 41 value 0.807374
## iter 42 value 0.806551
## iter 43 value 0.804976
## iter 44 value 0.803021
## iter 45 value 0.801629
## iter 46 value 0.801293
## iter 47 value 0.801083
## iter 48 value 0.801062
## iter 49 value 0.801057
## iter 50 value 0.801057
## iter 50 value 0.801057
## iter 50 value 0.801057
## final value 0.801057
## converged
## initial value 0.800808
## iter 2 value 0.797917
## iter 3 value 0.797560
## iter 4 value 0.797310
## iter 5 value 0.796981
## iter 6 value 0.796822
## iter 7 value 0.796772
## iter 8 value 0.796749
## iter 9 value 0.796746
## iter 10 value 0.796746
## iter 11 value 0.796745
## iter 12 value 0.796745
## iter 13 value 0.796744
## iter 14 value 0.796742
## iter 15 value 0.796739
## iter 16 value 0.796733
## iter 17 value 0.796724
## iter 18 value 0.796703
## iter 19 value 0.796700
## iter 20 value 0.796697
## iter 21 value 0.796695
## iter 22 value 0.796694
## iter 23 value 0.796694
## iter 24 value 0.796694
## iter 24 value 0.796694
## iter 24 value 0.796694
## final value 0.796694
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 sar1 sar2
## -0.0638 0.0922 0.3647 -0.0568 0.0504 -0.7192 0.7350 -0.5457
## s.e. 0.4768 0.3755 0.2256 0.0499 0.0501 0.4753 0.0664 0.0644
## sma1 sma2
## -1.2922 0.5803
## s.e. 0.0700 0.0604
##
## sigma^2 estimated as 4.752: log likelihood = -1021.41, aic = 2064.81
##
## $degrees_of_freedom
## [1] 451
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.0638 0.4768 -0.1338 0.8936
## ar2 0.0922 0.3755 0.2457 0.8060
## ar3 0.3647 0.2256 1.6165 0.1067
## ar4 -0.0568 0.0499 -1.1375 0.2560
## ar5 0.0504 0.0501 1.0067 0.3146
## ma1 -0.7192 0.4753 -1.5132 0.1309
## sar1 0.7350 0.0664 11.0656 0.0000
## sar2 -0.5457 0.0644 -8.4754 0.0000
## sma1 -1.2922 0.0700 -18.4594 0.0000
## sma2 0.5803 0.0604 9.6128 0.0000
##
## $AIC
## [1] 4.374604
##
## $AICc
## [1] 4.375615
##
## $BIC
## [1] 4.470933fit ARIMA by hand
plot(diff(df22, 3)) # d=1
acf(diff(df22, 3), lag.max=200) # MA:q= 0
pacf(diff(df22, 3), lag.max=200) # AR:p= 1, 2
# try ARIMA(1,1,0)
fit2.1 <- Arima(df22, order=c(1, 1, 0), seasonal=c(0,3,0), method="ML",include.drift = TRUE)## Warning in Arima(df22, order = c(1, 1, 0), seasonal = c(0, 3, 0), method =
## "ML", : No drift term fitted as the order of difference is 2 or more.summary(fit2.1)## Series: df22
## ARIMA(1,1,0)(0,3,0)[12]
##
## Coefficients:
## ar1
## -0.5182
## s.e. 0.0409
##
## sigma^2 estimated as 69.66: log likelihood=-1553.69
## AIC=3111.38 AICc=3111.4 BIC=3119.54
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.03230883 8.004988 5.838966 0.05872132 5.806634 2.023926
## ACF1
## Training set -0.2282839checkresiduals(fit2.1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(0,3,0)[12]
## Q* = 394.41, df = 23, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 24# AIC=3111.38 AICc=3111.4 BIC=3119.54
fit2.2 <- Arima(df22, order=c(2, 1, 0),seasonal=c(0,3,0), method="ML",include.drift = TRUE)## Warning in Arima(df22, order = c(2, 1, 0), seasonal = c(0, 3, 0), method =
## "ML", : No drift term fitted as the order of difference is 2 or more.summary(fit2.2)## Series: df22
## ARIMA(2,1,0)(0,3,0)[12]
##
## Coefficients:
## ar1 ar2
## -0.7483 -0.4426
## s.e. 0.0430 0.0430
##
## sigma^2 estimated as 56.21: log likelihood=-1506.52
## AIC=3019.05 AICc=3019.1 BIC=3031.29
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02702874 7.18248 5.234064 0.05320192 5.223328 1.814252
## ACF1
## Training set 0.02938732# AIC=3019.05 AICc=3019.1 BIC=3031.29
checkresiduals(fit2.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0)(0,3,0)[12]
## Q* = 269.44, df = 22, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 24# comparing the diagnostic plot, the ARIMA(5,1,1)(2,1,2)[12] is better in residual acfUse ARIMA model to predict future index
pred22=forecast(fit22,h=48)
pred22 # The last 2 columns are 95% prediction intervals lower bound and upper bound ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2019 97.15872 94.33405 99.98339 92.83876 101.47868
## Aug 2019 92.65677 89.76634 95.54720 88.23624 97.07730
## Sep 2019 97.44523 94.38186 100.50860 92.76021 102.13025
## Oct 2019 102.04632 98.48557 105.60707 96.60062 107.49202
## Nov 2019 101.66363 98.00648 105.32078 96.07050 107.25676
## Dec 2019 96.12563 92.21497 100.03629 90.14480 102.10646
## Jan 2020 96.89099 92.73152 101.05046 90.52963 103.25235
## Feb 2020 96.56224 92.26413 100.86035 89.98885 103.13564
## Mar 2020 106.70315 102.18151 111.22480 99.78789 113.61841
## Apr 2020 96.32467 91.62434 101.02500 89.13614 103.51321
## May 2020 97.71657 92.85921 102.57394 90.28788 105.14527
## Jun 2020 98.12829 93.08535 103.17123 90.41577 105.84080
## Jul 2020 96.24240 90.60779 101.87702 87.62501 104.85980
## Aug 2020 91.05382 85.20889 96.89875 82.11478 99.99287
## Sep 2020 99.35485 93.24877 105.46094 90.01640 108.69330
## Oct 2020 100.70513 94.25563 107.15463 90.84147 110.56878
## Nov 2020 101.64047 94.97898 108.30196 91.45260 111.82834
## Dec 2020 97.51976 90.58919 104.45032 86.92038 108.11913
## Jan 2021 96.47777 89.28843 103.66710 85.48263 107.47290
## Feb 2021 97.05218 89.64624 104.45811 85.72578 108.37858
## Mar 2021 107.13560 99.48614 114.78505 95.43676 118.83443
## Apr 2021 96.68426 88.81219 104.55632 84.64496 108.72355
## May 2021 97.77658 89.69261 105.86054 85.41322 110.13994
## Jun 2021 100.31640 92.01391 108.61888 87.61884 113.01395
## Jul 2021 95.93933 87.38900 104.48967 82.86272 109.01594
## Aug 2021 92.36534 83.60558 101.12509 78.96845 105.76223
## Sep 2021 101.05943 92.08641 110.03246 87.33638 114.78249
## Oct 2021 101.43646 92.24749 110.62543 87.38315 115.48977
## Nov 2021 103.85441 94.46598 113.24285 89.49605 118.21278
## Dec 2021 98.45137 88.86007 108.04267 83.78274 113.12000
## Jan 2022 98.19067 88.40147 107.97987 83.21938 113.16196
## Feb 2022 98.30338 88.32342 108.28333 83.04036 113.56640
## Mar 2022 108.46987 98.29905 118.64069 92.91495 124.02479
## Apr 2022 97.18001 86.82355 107.53647 81.34118 113.01884
## May 2022 99.06753 88.52899 109.60606 82.95023 115.18482
## Jun 2022 101.41191 90.69311 112.13071 85.01892 117.80490
## Jul 2022 96.39425 85.44720 107.34130 79.65218 113.13631
## Aug 2022 94.33364 83.20311 105.46418 77.31096 111.35633
## Sep 2022 101.42528 90.10700 112.74355 84.11547 118.73508
## Oct 2022 102.86239 91.34590 114.37888 85.24944 120.47534
## Nov 2022 105.63436 93.93951 117.32921 87.74863 123.52009
## Dec 2022 98.53044 86.65072 110.41015 80.36199 116.69888
## Jan 2023 99.82446 87.76232 111.88660 81.37700 118.27191
## Feb 2023 99.10189 86.86534 111.33844 80.38771 117.81607
## Mar 2023 109.36698 96.95338 121.78058 90.38202 128.35195
## Apr 2023 97.49608 84.90984 110.08231 78.24709 116.74506
## May 2023 100.13238 87.37672 112.88805 80.62427 119.64049
## Jun 2023 101.17348 88.24863 114.09834 81.40662 120.94034autoplot(pred22) #plot with confidence interval 
(4) ARCH/GARCH model
# fit ARIMA model without seasonal term
fit3=auto.arima(df22, seasonal = FALSE, trace = TRUE) # ARIMA(2,1,2)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : 2911.535
## ARIMA(0,1,0) with drift : 3129.63
## ARIMA(1,1,0) with drift : 3081.236
## ARIMA(0,1,1) with drift : 2939.252
## ARIMA(0,1,0) : 3127.617
## ARIMA(1,1,2) with drift : 2934.21
## ARIMA(2,1,1) with drift : 2926.379
## ARIMA(3,1,2) with drift : Inf
## ARIMA(2,1,3) with drift : 2832.569
## ARIMA(1,1,3) with drift : 2934.936
## ARIMA(3,1,3) with drift : Inf
## ARIMA(2,1,4) with drift : Inf
## ARIMA(1,1,4) with drift : Inf
## ARIMA(3,1,4) with drift : Inf
## ARIMA(2,1,3) : 2831.025
## ARIMA(1,1,3) : 2933.27
## ARIMA(2,1,2) : 2909.949
## ARIMA(3,1,3) : Inf
## ARIMA(2,1,4) : Inf
## ARIMA(1,1,2) : 2932.669
## ARIMA(1,1,4) : Inf
## ARIMA(3,1,2) : Inf
## ARIMA(3,1,4) : Inf
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,1,3) : Inf
## ARIMA(2,1,3) with drift : Inf
## ARIMA(2,1,2) : 2907.698
##
## Best model: ARIMA(2,1,2)df22.diff=diff(df22)
plot(df22.diff, main=" diff data") 
df22.log=log(df22)
plot(df22.log, main=" log data")
di22.difflog=diff(log(df22))
plot(di22.difflog, main=" difflog data")
acf(di22.difflog, lag.max = 100) # from pacf plot, we could observe the seasonality
pacf(di22.difflog, lag.max = 100)# from pacf plot, AR(3) or AR(4)
acf2(df22^2) # accoring to pacf plot, ARCH(4) 
## ACF PACF
## [1,] 0.78 0.78
## [2,] 0.71 0.26
## [3,] 0.72 0.30
## [4,] 0.76 0.32
## [5,] 0.73 0.12
## [6,] 0.75 0.24
## [7,] 0.71 -0.03
## [8,] 0.75 0.25
## [9,] 0.68 -0.17
## [10,] 0.65 -0.08
## [11,] 0.75 0.35
## [12,] 0.90 0.64
## [13,] 0.71 -0.45
## [14,] 0.65 -0.32
## [15,] 0.62 -0.41
## [16,] 0.68 0.05
## [17,] 0.66 0.02
## [18,] 0.65 0.08
## [19,] 0.62 0.09
## [20,] 0.68 0.09
## [21,] 0.58 -0.04
## [22,] 0.57 0.09
## [23,] 0.67 0.07
## [24,] 0.78 0.06
## [25,] 0.63 -0.01
## [26,] 0.57 -0.18
## [27,] 0.52 -0.13
## [28,] 0.60 0.00
## [29,] 0.57 -0.05
## [30,] 0.56 0.01
## [31,] 0.55 0.11
## [32,] 0.58 0.03
## [33,] 0.48 0.01
## [34,] 0.50 0.06
## [35,] 0.57 -0.16
## [36,] 0.68 0.09
## [37,] 0.56 0.08
## [38,] 0.47 -0.12
## [39,] 0.43 -0.10
## [40,] 0.53 0.05
## [41,] 0.47 0.00
## [42,] 0.47 -0.03
## [43,] 0.47 0.08
## [44,] 0.48 -0.02
## [45,] 0.40 -0.05
## [46,] 0.43 0.13
## [47,] 0.46 -0.14
## [48,] 0.60 0.08# fit ARIMA(2,1,2) model, then plot the residuals ACF and PACF
arima212.2=arima(df22, order=c(2,1,2))
residual_arima212.2=arima212.2$residuals
acf2(residual_arima212.2^2)
## ACF PACF
## [1,] -0.04 -0.04
## [2,] -0.20 -0.20
## [3,] -0.14 -0.16
## [4,] 0.07 0.01
## [5,] 0.09 0.04
## [6,] -0.09 -0.10
## [7,] 0.15 0.19
## [8,] 0.09 0.10
## [9,] -0.13 -0.09
## [10,] -0.14 -0.08
## [11,] -0.11 -0.16
## [12,] 0.52 0.45
## [13,] 0.01 0.02
## [14,] -0.17 -0.04
## [15,] -0.10 0.00
## [16,] 0.04 -0.01
## [17,] 0.07 0.02
## [18,] -0.12 -0.04
## [19,] 0.20 0.13
## [20,] 0.06 -0.05
## [21,] -0.10 0.02
## [22,] -0.18 -0.10
## [23,] -0.09 -0.02
## [24,] 0.39 0.13
## [25,] 0.03 -0.02
## [26,] -0.19 -0.11
## [27,] -0.05 0.03
## [28,] -0.01 -0.06
## [29,] 0.12 0.09
## [30,] -0.11 0.04
## [31,] 0.20 0.04
## [32,] 0.05 -0.01
## [33,] -0.06 0.05
## [34,] -0.19 -0.07
## [35,] -0.04 0.06
## [36,] 0.29 0.03
## [37,] 0.00 -0.09
## [38,] -0.18 -0.05
## [39,] -0.07 -0.07
## [40,] 0.02 0.03
## [41,] 0.14 0.10
## [42,] -0.12 -0.02
## [43,] 0.15 -0.01
## [44,] 0.04 0.00
## [45,] -0.10 -0.07
## [46,] -0.19 -0.06
## [47,] -0.05 0.01
## [48,] 0.42 0.22par(mfrow=c(2,2))
acf(residual_arima212.2, main="ACF of residuals (UK)", lag.max=200)
pacf(residual_arima212.2, main="PACF of residuals (UK)", lag.max=200)
acf(residual_arima212.2^2, main="ACF of squared residuals (UK)", lag.max=200)
pacf(residual_arima212.2^2, main="PACF of squared residuals (UK)", lag.max=200)
garch_uk1=garch(residual_arima212,order=c(0,7),trace=F)
summary(garch_uk1)##
## Call:
## garch(x = residual_arima212, order = c(0, 7), trace = F)
##
## Model:
## GARCH(0,7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.75116 -0.57509 0.05721 0.79602 2.90765
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 1.225e+00 2.166e-01 5.655 1.56e-08 ***
## a1 5.160e-02 6.214e-02 0.830 0.406
## a2 6.612e-02 6.110e-02 1.082 0.279
## a3 5.267e-02 5.364e-02 0.982 0.326
## a4 1.972e-02 6.305e-02 0.313 0.754
## a5 5.674e-02 6.206e-02 0.914 0.361
## a6 4.078e-15 5.792e-02 0.000 1.000
## a7 7.533e-02 5.892e-02 1.278 0.201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 7.5789, df = 2, p-value = 0.02261
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.0044309, df = 1, p-value = 0.9469AIC(garch_uk1) # the AIC is smaller, so ARCH(0,7) model is better ## [1] 1641.352garch_uk2=garch(residual_arima212,order=c(4,7),trace=F) # (q,p)
summary(garch_uk2)##
## Call:
## garch(x = residual_arima212, order = c(4, 7), trace = F)
##
## Model:
## GARCH(4,7)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.71673 -0.57249 0.05596 0.78096 2.81972
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 8.502e-01 9.173e-01 0.927 0.354
## a1 4.923e-02 6.722e-02 0.732 0.464
## a2 6.343e-02 7.307e-02 0.868 0.385
## a3 4.842e-02 8.977e-02 0.539 0.590
## a4 1.470e-02 8.829e-02 0.166 0.868
## a5 5.396e-02 7.459e-02 0.723 0.469
## a6 3.423e-15 8.823e-02 0.000 1.000
## a7 8.089e-02 8.127e-02 0.995 0.320
## b1 5.343e-02 8.972e-01 0.060 0.953
## b2 5.570e-02 8.293e-01 0.067 0.946
## b3 5.947e-02 5.870e-01 0.101 0.919
## b4 6.267e-02 5.587e-01 0.112 0.911
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 7.2065, df = 2, p-value = 0.02723
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.019276, df = 1, p-value = 0.8896AIC(garch_uk2) ## [1] 1649.505Canada
Research Indicator: Economic Activity, Industrial Production, Index
Indicator Code: AIP_IX
Country: Canada
df1_ca=subset(ifs, (Indicator_Code=='AIP_IX') & (Country_Name=='Canada'))
head(df1_ca)## Country_Name Country_Code
## 178902 Canada 156
## 178903 Canada 156
## 178904 Canada 156
## 178905 Canada 156
## 178906 Canada 156
## 178907 Canada 156
## Indicator_Name Indicator_Code
## 178902 Economic Activity, Industrial Production, Index AIP_IX
## 178903 Economic Activity, Industrial Production, Index AIP_IX
## 178904 Economic Activity, Industrial Production, Index AIP_IX
## 178905 Economic Activity, Industrial Production, Index AIP_IX
## 178906 Economic Activity, Industrial Production, Index AIP_IX
## 178907 Economic Activity, Industrial Production, Index AIP_IX
## Time_Period Value
## 178902 1997M1 NA
## 178903 2007M5 108.47337
## 178904 2007M8 108.61256
## 178905 2008M4 105.83595
## 178906 2009M1 99.02647
## 178907 2009M4 92.17060change the format of date
df1_ca$Time_Period <- gsub("M", "/", df1_ca$Time_Period)
df1_ca['date'] = '/01'
df1_ca$Time_Period <- paste(df1_ca$Time_Period, df1_ca$date, sep="")
df1_ca$Time_Period = as.Date(df1_ca$Time_Period, "%Y/%m/%d")
df1_ca = df1_ca[order(df1_ca$Time_Period),]
df1_ca = subset(df1_ca, select = -c(date))
rownames(df1_ca)=NULL
df1_ca=na.omit(df1_ca)
head(df1_ca)## Country_Name Country_Code
## 2 Canada 156
## 3 Canada 156
## 4 Canada 156
## 5 Canada 156
## 6 Canada 156
## 7 Canada 156
## Indicator_Name Indicator_Code
## 2 Economic Activity, Industrial Production, Index AIP_IX
## 3 Economic Activity, Industrial Production, Index AIP_IX
## 4 Economic Activity, Industrial Production, Index AIP_IX
## 5 Economic Activity, Industrial Production, Index AIP_IX
## 6 Economic Activity, Industrial Production, Index AIP_IX
## 7 Economic Activity, Industrial Production, Index AIP_IX
## Time_Period Value
## 2 2000-01-01 106.4326
## 3 2000-02-01 112.6132
## 4 2000-03-01 114.5484
## 5 2000-04-01 110.3226
## 6 2000-05-01 112.6478
## 7 2000-06-01 115.5850Only keep value and Time_Period
# Research Indicator: Economic Activity, Industrial Production, Index
# Indicator Code: AIP_IX
# Country: Canada
# the original data
df33=ts(df1_ca$Value, frequency = 12, start = c(2000, 1), end=c(2018, 12))
plot(df33, main="Canada: Industrial Production Index")
(1) decompose Analysis
plot(decompose(df33))
(2) Spectral Analysis
sp=spectrum(df33, log='no')
# find the cycle with predominant frequecny
sp$fre[which.max(sp$spec)] # 0.05## [1] 0.05# find the cycle of predominant frequecny
1/sp$fre[which.max(sp$spec)] # 20 month## [1] 20max(sp$spec) #estimate of spectral density at predominant period which is close to 115.1996.## [1] 115.1996U = qchisq(.025,2)
L = qchisq(.975,2)
2*max(sp$spec)/L #lower bound for spectral density## [1] 31.228882*max(sp$spec)/U #upper bound for spectral density## [1] 4550.14# The 95% confidence interval is [31.22888, 4550.14](3) ARIMA model
fit33 <- auto.arima(df33, trace = TRUE) # Best model: ARIMA(2,1,2)(0,1,1)[12]##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2)(1,1,1)[12] : 812.4419
## ARIMA(0,1,0)(0,1,0)[12] : 859.5941
## ARIMA(1,1,0)(1,1,0)[12] : 820.3644
## ARIMA(0,1,1)(0,1,1)[12] : 815.9909
## ARIMA(2,1,2)(0,1,1)[12] : 805.7274
## ARIMA(2,1,2)(0,1,0)[12] : 855.5499
## ARIMA(2,1,2)(0,1,2)[12] : 807.6211
## ARIMA(2,1,2)(1,1,0)[12] : 817.1942
## ARIMA(2,1,2)(1,1,2)[12] : 809.6784
## ARIMA(1,1,2)(0,1,1)[12] : 817.9261
## ARIMA(2,1,1)(0,1,1)[12] : 813.2687
## ARIMA(3,1,2)(0,1,1)[12] : 810.3778
## ARIMA(2,1,3)(0,1,1)[12] : 811.8729
## ARIMA(1,1,1)(0,1,1)[12] : 818.3023
## ARIMA(1,1,3)(0,1,1)[12] : 813.2287
## ARIMA(3,1,1)(0,1,1)[12] : 809.6983
## ARIMA(3,1,3)(0,1,1)[12] : 811.5295
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(2,1,2)(0,1,1)[12] : 825.5495
##
## Best model: ARIMA(2,1,2)(0,1,1)[12]evaluate the model
sarima(df33, 2,1,2,0,1,1,12) # good fit ## initial value 0.611630
## iter 2 value 0.514917
## iter 3 value 0.502950
## iter 4 value 0.497509
## iter 5 value 0.494460
## iter 6 value 0.494060
## iter 7 value 0.493732
## iter 8 value 0.492755
## iter 9 value 0.492085
## iter 10 value 0.491685
## iter 11 value 0.491660
## iter 12 value 0.491657
## iter 13 value 0.491655
## iter 14 value 0.491633
## iter 15 value 0.491286
## iter 16 value 0.491085
## iter 17 value 0.490687
## iter 18 value 0.490617
## iter 19 value 0.490482
## iter 20 value 0.490400
## iter 21 value 0.490255
## iter 22 value 0.487365
## iter 23 value 0.485588
## iter 24 value 0.485003
## iter 25 value 0.484367
## iter 26 value 0.482405
## iter 27 value 0.481542
## iter 28 value 0.480885
## iter 29 value 0.475778
## iter 30 value 0.473318
## iter 31 value 0.472775
## iter 32 value 0.471360
## iter 33 value 0.470405
## iter 34 value 0.468233
## iter 35 value 0.467115
## iter 36 value 0.465970
## iter 37 value 0.465215
## iter 38 value 0.464134
## iter 39 value 0.463384
## iter 40 value 0.463122
## iter 41 value 0.463061
## iter 42 value 0.463048
## iter 43 value 0.463037
## iter 44 value 0.463032
## iter 45 value 0.463032
## iter 46 value 0.463032
## iter 46 value 0.463032
## iter 46 value 0.463032
## final value 0.463032
## converged
## initial value 0.475590
## iter 2 value 0.473109
## iter 3 value 0.472911
## iter 4 value 0.472845
## iter 5 value 0.472148
## iter 6 value 0.472138
## iter 7 value 0.472136
## iter 8 value 0.472135
## iter 9 value 0.472127
## iter 10 value 0.472115
## iter 11 value 0.472101
## iter 12 value 0.472098
## iter 13 value 0.472098
## iter 14 value 0.472098
## iter 15 value 0.472098
## iter 16 value 0.472098
## iter 17 value 0.472098
## iter 17 value 0.472098
## iter 17 value 0.472098
## final value 0.472098
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 ma2 sma1
## -1.1088 -0.8565 0.9736 0.7177 -0.4960
## s.e. 0.1201 0.1754 0.1805 0.2449 0.0568
##
## sigma^2 estimated as 2.526: log likelihood = -406.57, aic = 825.15
##
## $degrees_of_freedom
## [1] 210
##
## $ttable
## Estimate SE t.value p.value
## ar1 -1.1088 0.1201 -9.2346 0.0000
## ar2 -0.8565 0.1754 -4.8837 0.0000
## ma1 0.9736 0.1805 5.3937 0.0000
## ma2 0.7177 0.2449 2.9306 0.0038
## sma1 -0.4960 0.0568 -8.7308 0.0000
##
## $AIC
## [1] 3.651087
##
## $AICc
## [1] 3.652294
##
## $BIC
## [1] 3.740573fit ARIMA by hand
plot(diff(df33,4)) # d=1
acf(diff(df33,4), lag.max=200) # MA:q= 0
pacf(diff(df33,4), lag.max=200) # AR:p= 1, 2
# try ARIMA(1,1,0)
fit3.1 <- Arima(df33, order=c(1, 1, 0), seasonal=c(0,4,0), method="ML",include.drift = TRUE)## Warning in Arima(df33, order = c(1, 1, 0), seasonal = c(0, 4, 0), method =
## "ML", : No drift term fitted as the order of difference is 2 or more.summary(fit3.1)## Series: df33
## ARIMA(1,1,0)(0,4,0)[12]
##
## Coefficients:
## ar1
## -0.0911
## s.e. 0.0743
##
## sigma^2 estimated as 95.66: log likelihood=-697.4
## AIC=1398.8 AICc=1398.86 BIC=1405.17
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05204617 8.642069 5.868094 0.06058388 5.41006 1.546495
## ACF1
## Training set 0.01003916checkresiduals(fit3.1)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0)(0,4,0)[12]
## Q* = 164.83, df = 23, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 24# AIC=3111.38 AICc=3111.4 BIC=3119.54
fit3.2 <- Arima(df33, order=c(2, 1, 0),seasonal=c(0,4,0), method="ML",include.drift = TRUE)## Warning in Arima(df33, order = c(2, 1, 0), seasonal = c(0, 4, 0), method =
## "ML", : No drift term fitted as the order of difference is 2 or more.summary(fit3.2)## Series: df33
## ARIMA(2,1,0)(0,4,0)[12]
##
## Coefficients:
## ar1 ar2
## -0.0814 0.1074
## s.e. 0.0741 0.0740
##
## sigma^2 estimated as 95.06: log likelihood=-696.35
## AIC=1398.71 AICc=1398.84 BIC=1408.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04387278 8.590645 5.803412 0.05248819 5.349517 1.529449
## ACF1
## Training set -0.01498271# AIC=3019.05 AICc=3019.1 BIC=3031.29
checkresiduals(fit3.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0)(0,4,0)[12]
## Q* = 155.59, df = 22, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 24Use ARIMA model to predict future index
pred33=forecast(fit33,h=48)
pred33 # The last 2 columns are 95% prediction intervals lower bound and upper bound ## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2019 122.6092 120.5480 124.6704 119.45692 125.7615
## Feb 2019 126.6568 123.9318 129.3817 122.48926 130.8243
## Mar 2019 125.0077 121.7389 128.2765 120.00851 130.0068
## Apr 2019 119.5537 115.7118 123.3955 113.67801 125.4293
## May 2019 118.8311 114.6042 123.0579 112.36661 125.2955
## Jun 2019 121.7887 117.1692 126.4083 114.72371 128.8538
## Jul 2019 119.1687 114.1463 124.1912 111.48755 126.8499
## Aug 2019 122.8566 117.5329 128.1803 114.71470 130.9984
## Sep 2019 123.6803 118.0272 129.3334 115.03459 132.3259
## Oct 2019 124.2781 118.3057 130.2505 115.14408 133.4120
## Nov 2019 128.5006 122.2679 134.7332 118.96859 138.0326
## Dec 2019 122.2585 115.7366 128.7804 112.28409 132.2329
## Jan 2020 125.3760 118.2263 132.5257 114.44144 136.3105
## Feb 2020 129.2488 121.6050 136.8927 117.55853 140.9391
## Mar 2020 127.9106 119.7576 136.0635 115.44173 140.3794
## Apr 2020 122.2615 113.6139 130.9092 109.03613 135.4870
## May 2020 121.4890 112.4244 130.5535 107.62585 135.3520
## Jun 2020 124.6691 115.1639 134.1743 110.13210 139.2061
## Jul 2020 121.8452 111.9237 131.7667 106.67164 137.0188
## Aug 2020 125.5686 115.2757 135.8614 109.82699 141.3101
## Sep 2020 126.5275 115.8422 137.2129 110.18572 142.8693
## Oct 2020 126.9449 115.8951 137.9948 110.04570 143.8442
## Nov 2020 131.2516 119.8612 142.6420 113.83151 148.6717
## Dec 2020 125.0707 113.3251 136.8164 107.10730 143.0342
## Jan 2021 128.0483 115.6924 140.4042 109.15152 146.9450
## Feb 2021 132.0238 119.1301 144.9176 112.30449 151.7432
## Mar 2021 130.6916 117.2585 144.1246 110.14750 151.2356
## Apr 2021 124.9479 110.9865 138.9094 103.59574 146.3001
## May 2021 124.2751 109.8335 138.7167 102.18855 146.3617
## Jun 2021 127.4256 112.4962 142.3551 104.59300 150.2583
## Jul 2021 124.5491 109.1489 139.9494 100.99654 148.1018
## Aug 2021 128.3562 112.5161 144.1963 104.13085 152.5815
## Sep 2021 129.2674 112.9797 145.5551 104.35754 154.1773
## Oct 2021 129.6661 112.9502 146.3819 104.10139 155.2308
## Nov 2021 134.0344 116.9090 151.1598 107.84340 160.2255
## Dec 2021 127.8012 110.2612 145.3413 100.97604 154.6264
## Jan 2022 130.7840 112.6029 148.9651 102.97838 158.5895
## Feb 2022 134.7986 116.0358 153.5614 106.10334 163.4939
## Mar 2022 133.4185 114.0770 152.7600 103.83824 162.9988
## Apr 2022 127.6944 107.7764 147.6124 97.23247 158.1563
## May 2022 127.0409 106.5901 147.4917 95.76409 158.3177
## Jun 2022 130.1533 109.1654 151.1412 98.05512 162.2515
## Jul 2022 127.3025 105.7879 148.8172 94.39878 160.2063
## Aug 2022 131.1137 109.1028 153.1246 97.45093 164.7765
## Sep 2022 131.9983 109.4852 154.5115 97.56747 166.4292
## Oct 2022 132.4230 109.4223 155.4236 97.24645 167.5995
## Nov 2022 136.7853 113.3170 160.2536 100.89370 172.6769
## Dec 2022 130.5365 106.5960 154.4770 93.92272 167.1503autoplot(pred33) #plot with confidence interval 
(4) ARCH/GARCH model
# fit ARIMA model without seasonal term
fit4=auto.arima(df33, seasonal = FALSE, trace = TRUE) # ARIMA(1,1,2)##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,1,2) with drift : 1265.845
## ARIMA(0,1,0) with drift : 1338.367
## ARIMA(1,1,0) with drift : 1318.347
## ARIMA(0,1,1) with drift : 1287.103
## ARIMA(0,1,0) : 1336.366
## ARIMA(1,1,2) with drift : 1262.375
## ARIMA(0,1,2) with drift : 1279.524
## ARIMA(1,1,1) with drift : 1280.805
## ARIMA(1,1,3) with drift : 1264.084
## ARIMA(0,1,3) with drift : 1271.954
## ARIMA(2,1,1) with drift : 1279.615
## ARIMA(2,1,3) with drift : 1267.933
## ARIMA(1,1,2) : 1260.591
## ARIMA(0,1,2) : 1277.917
## ARIMA(1,1,1) : 1278.922
## ARIMA(2,1,2) : 1263.951
## ARIMA(1,1,3) : 1262.254
## ARIMA(0,1,1) : 1285.394
## ARIMA(0,1,3) : 1270.212
## ARIMA(2,1,1) : 1277.722
## ARIMA(2,1,3) : 1266.006
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(1,1,2) : 1265.904
##
## Best model: ARIMA(1,1,2)df33.diff=diff(df33)
plot(df33.diff, main=" diff data") 
df33.log=log(df33)
plot(df33.log, main=" log data")
di33.difflog=diff(log(df33))
plot(di33.difflog, main=" difflog data")
acf(di33.difflog, lag.max = 100) # from pacf plot, we could observe the seasonality
pacf(di33.difflog, lag.max = 100)# from pacf plot, AR(2) or AR(4)
acf2(df33^2) # accoring to pacf plot, ARCH(2) 
## ACF PACF
## [1,] 0.76 0.76
## [2,] 0.65 0.18
## [3,] 0.66 0.30
## [4,] 0.62 0.04
## [5,] 0.63 0.22
## [6,] 0.58 -0.07
## [7,] 0.57 0.15
## [8,] 0.52 -0.14
## [9,] 0.51 0.14
## [10,] 0.45 -0.26
## [11,] 0.49 0.41
## [12,] 0.66 0.32
## [13,] 0.44 -0.65
## [14,] 0.33 -0.07
## [15,] 0.34 -0.01
## [16,] 0.30 0.11
## [17,] 0.32 -0.08
## [18,] 0.27 0.04
## [19,] 0.27 0.07
## [20,] 0.23 -0.03
## [21,] 0.23 0.12
## [22,] 0.19 0.12
## [23,] 0.24 0.02
## [24,] 0.39 -0.07
## [25,] 0.19 -0.13
## [26,] 0.10 -0.08
## [27,] 0.11 -0.03
## [28,] 0.10 0.03
## [29,] 0.12 -0.07
## [30,] 0.08 0.13
## [31,] 0.10 0.00
## [32,] 0.08 0.10
## [33,] 0.09 -0.01
## [34,] 0.06 -0.02
## [35,] 0.10 -0.01
## [36,] 0.25 0.07
## [37,] 0.08 -0.15
## [38,] -0.01 -0.03
## [39,] 0.01 -0.03
## [40,] -0.01 -0.09
## [41,] 0.01 -0.06
## [42,] -0.03 0.00
## [43,] -0.03 0.02
## [44,] -0.05 -0.05
## [45,] -0.03 0.06
## [46,] -0.08 0.00
## [47,] -0.05 -0.04
## [48,] 0.08 -0.06# fit ARIMA(1,1,2) model, then plot the residuals ACF and PACF
arima112=arima(df33, order=c(1,1,2))
residual_arima112=arima112$residuals
acf2(residual_arima112^2)
## ACF PACF
## [1,] -0.05 -0.05
## [2,] -0.06 -0.06
## [3,] -0.03 -0.04
## [4,] -0.09 -0.10
## [5,] 0.28 0.27
## [6,] -0.17 -0.17
## [7,] 0.24 0.29
## [8,] -0.08 -0.14
## [9,] -0.03 0.11
## [10,] -0.07 -0.27
## [11,] -0.06 0.20
## [12,] 0.68 0.55
## [13,] -0.06 0.02
## [14,] -0.07 -0.11
## [15,] -0.05 0.03
## [16,] -0.08 -0.05
## [17,] 0.24 0.01
## [18,] -0.14 -0.02
## [19,] 0.20 0.08
## [20,] -0.09 -0.09
## [21,] -0.02 0.07
## [22,] -0.02 0.03
## [23,] -0.08 0.01
## [24,] 0.50 -0.03
## [25,] -0.09 0.00
## [26,] -0.07 -0.07
## [27,] -0.02 0.10
## [28,] -0.04 0.01
## [29,] 0.17 -0.04
## [30,] -0.12 -0.04
## [31,] 0.17 0.06
## [32,] -0.05 0.02
## [33,] -0.01 0.03
## [34,] -0.02 -0.06
## [35,] -0.07 0.04
## [36,] 0.47 0.19
## [37,] -0.06 0.07
## [38,] -0.07 -0.04
## [39,] -0.03 -0.06
## [40,] -0.04 -0.04
## [41,] 0.13 -0.04
## [42,] -0.10 0.04
## [43,] 0.16 0.05
## [44,] -0.02 0.02
## [45,] 0.01 0.02
## [46,] -0.02 0.06
## [47,] -0.07 -0.01
## [48,] 0.38 -0.03par(mfrow=c(2,2))
acf(residual_arima112, main="ACF of residuals (Canada)", lag.max=100)
pacf(residual_arima112, main="PACF of residuals (Canada)", lag.max=100)
acf(residual_arima112^2, main="ACF of squared residuals (Canada)", lag.max=100)
pacf(residual_arima112^2, main="PACF of squared residuals (Canada)", lag.max=100)
garch_ca1=garch(residual_arima212,order=c(6,1),trace=F)
summary(garch_ca1)##
## Call:
## garch(x = residual_arima212, order = c(6, 1), trace = F)
##
## Model:
## GARCH(6,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.72773 -0.55433 0.05367 0.78023 2.64293
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 1.204e+00 2.043e+00 0.590 0.555
## a1 6.910e-02 7.946e-02 0.870 0.385
## b1 7.853e-02 4.716e-01 0.167 0.868
## b2 6.655e-02 4.747e-01 0.140 0.889
## b3 1.636e-14 8.518e-01 0.000 1.000
## b4 6.499e-02 6.449e-01 0.101 0.920
## b5 6.557e-03 7.233e-01 0.009 0.993
## b6 8.502e-02 7.043e-01 0.121 0.904
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 9.3336, df = 2, p-value = 0.009402
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.47876, df = 1, p-value = 0.489AIC(garch_ca1) ## [1] 1656.202garch_ca2=garch(residual_arima112,order=c(6,2),trace=F) # (q,p) # AIC is smaller
summary(garch_ca2)##
## Call:
## garch(x = residual_arima112, order = c(6, 2), trace = F)
##
## Model:
## GARCH(6,2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5582 -0.6079 0.2180 0.7475 1.7741
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 8.902e+00 1.150e+01 0.774 0.439
## a1 8.104e-15 9.045e-02 0.000 1.000
## a2 4.910e-03 9.147e-02 0.054 0.957
## b1 5.858e-02 2.163e+01 0.003 0.998
## b2 5.763e-02 1.612e+01 0.004 0.997
## b3 6.608e-02 6.839e+00 0.010 0.992
## b4 6.506e-02 1.920e+01 0.003 0.997
## b5 6.710e-02 9.518e+00 0.007 0.994
## b6 6.719e-02 1.964e+01 0.003 0.997
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 26.185, df = 2, p-value = 2.061e-06
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.54673, df = 1, p-value = 0.4597AIC(garch_ca2) # the AIC is smaller, so GARCH(6,2) model is better ## [1] 1248.765